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ABSTRACT 

In  this  thesis,  the  use  of  the  Wentzel-Kramers-Brillouin  (WKB)  Theory  to  obtain  the 
solution  to  the  Helmholtz  Equation  governing  the  acoustic  normal  modes  is  examined. 
Specifically,  uniformly  valid  WKB  solutions  for  four  classes  of  acoustic  normal  modes 
in  the  ocean  are  derived  and  the  accuracy  of  the  WKB  approximation  is  tested  against 
some  exact  solutions.  It  is  found  that  this  inherently  high  frequency  technique  has  an 
appreciable  accuracy  even  at  a  frequency  of  1  Hz.  A  product  of  this  thesis  is  a  computer 
program  that  solves  for  the  WKB  modes  for  an  arbitrary  sound  speed  profile. 
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I.     INTRODUCTION 

A.  BACKGROUND 

To  model  sound  propagation  in  the  ocean  we  can  use  several  theories,  namely  Ray 
Theory,  Parabolic  Equation  Approximation  and  Normal  Modes.  The  Ray  Theory  sol- 
ution is  an  asymptotic  geometric-optics  solution,  obtainable  using  simple  ray  tracing 
techniques.  It  provides  a  simple  physical  description  on  how  sound  is  transmitted 
underwater.  However,  it  neglects  sound  diffraction  and  thus  needs  corrections  near 
caustics  and  turning  points.  Such  corrections  can  be  very  complicated  mathematically. 
The  Parabolic  Equation  method  is  less  physical  but  is  a  full-wave  solution.  An  asset  of 
Parabolic  Equation  is  its  capability  to  handle  variable  bottom  bathymetry  well.  Its  lim- 
itation is  that  it  does  not  accurately  model  sound  energy  propagating  in  steep  angles. 
Normal  Mode  Theory  gives  a  physical  full-wave  solution  to  the  wave  equation.  It  has 
some  computational  difficulties  in  handling  range  varying  sound  speed  fields. 

A  computational  difficulty  associated  with  normal  mode  models  applied  to  a  range- 
dependent  ocean  is  that  the  normal  modes  associated  with  a  great  number  of  profiles 
must  be  computed  and  the  results  stored.  Large  storage  and  processing  time  are  re- 
quired if  straight-forward  numerical  methods  such  as  finite  differences  are  used  to  com- 
pute the  normal  modes  (Chiu  and  Ehret,  1990).  The  use  of  finite  differences  methods 
requires  the  discretized  mode  functions  and  their  eigenvalues  at  a  great  number  of  points 
to  be  stored.  Moreover,  at  high  frequencies  these  methods  require  the  computation  of 
eigenvectors  and  eigenvalues  of  large  matrices,  which  can  result  in  significant  increases 
of  processing  time  and  numerical  noise  in  the  solution.  In  this  thesis  we  examine  an 
asymptotic  expansion  method  that  has  the  potential  to  overcome  this  difficulty.  The 
method  is  called  Wentzel-Kramers-Brillouin  (WKB)  theory. 

B.  THE  NORMAL  MODE  APPROACH 

The  wave  equation  governing  the  sound  pressure  p  in  the  ocean  is 

or 

where  c  =  c{z\  r,  6)  is  the  sound  speed  and  z  is  the  vertical  coordinate,  r  the  range  and  6 
the  azimuthal  angle.    In  the  three-dimensional  coupled  mode  model  of  Chiu  and  Ehret 


(1990),  the  pressure  is  expressed  as  a  linear  combination  of  the  local  normal  modes  Zn 
such  that 


p(r,d,z)  =  YJZn(z;r,d)Pn(r,e)T(t) 


For  a  harmonic  frequency  time  dependence 

T(t)  =  eia}1 

where  w  is  the  acoustic  angular  frequency,  the  local  modes  at  each  horizontal  location 
(r,  6)  are  required  to  satisfy 

<?2Z 


oz         L  J 


and  the  appropriate  boundary  conditions.  This  equation  is  usually  known  as  the 
Helmholtz  Equation.  The  constant  k„  is  the  horizontal  component  of  the  wavenumber 
vector  whose  magnitude  is  given  by 

CO 

k  =  —    . 

In  general,  there  are  many  possible  values  of  xn  (eigenvalues)  satisfying  the  Helmholtz 
Equation.  For  each  k„  there  is  an  associated  mode  Zn  (eigenfunction).  In  a  range- 
dependent  sound  speed  field  this  eigenvalue-eigenfunction  problem  must  be  solved  at  a 
great  number  of  horizontal  grid  points. 

C.     THESIS  OBJECTIVES  AND  OUTLINE 

The  main  objective  of  this  thesis  is  to  examine  the  use  of  the  WKB  theory  to  solve 
the  Helmholtz  Equation  governing  the  acoustic  normal  modes.  This  includes: 

1.  the  development  of  the  various  WKB  formulae  for  the  four  classes  of  normal  modes 
which  can  exist  in  single  duct; channel  environments, 

2.  the  parameterization  of  each  class  of  normal  modes  using  a  minimum  number  of 
parameters, 

3.  the  quantification  of  errors  in  the  WKB  solution  through  comparisons  to  three 
exact  solutions. 

A  product  coming  out  of  this  thesis  is  a  computer  program  solving  for  the  mode 
parameters  for  an  arbitrary  sound  speed  profile.  The  incorporation  of  this  code  in  the 


three-dimensional  coupled  mode  model  of  Chiu  and  Ehret  (1990)  is  expected  to  result  in 
significant  savings  of  processing  time  and  computer  storage. 

In  Chapter  II,  the  WKB  formulae  are  developed.  The  four  types  of  normal  modes 
are  discussed.  In  Chapter  III,  the  normal  modes  and  their  eigenvalues  for  three  abstract 
sound  speed  profiles  are  solved  exactly.  The  WKB  results  are  then  compared  to  the  exact 
solutions  for  an  error  analysis.  Conclusions  are  given  in  Chapter  IV. 


II.     WKB  NORMAL  MODES 

A.     FIRST  ORDER  WKB  APPROXIMATION 
1.     The  Mathematical  Problem 

The  WKB  approximation  to  the  solution  of  a  differential  equation  is  an 
asymptotic  expansion  method  applicable  to  a  large  range  of  equations.  It  is  named  after 
Wentzel,  Kramers  and  Brillouin  who  used  it  separately  but  at  about  the  same  time  in 
1926.  However  the  principle  of  this  technique  was  developed  by  Liouville  and  Green  in 
1837.  It  was  also  used  by  Rayleigh  (1912),  Gans  (1915)  and  Jeffreys  (1924),  among  oth- 
ers. The  WKB  method  is  also  known  as  the  Liouville-Green  or  WKBJ  approximation 
(the  letter  J  is  used  to  honour  Jeffreys'  contribution). 

To  compute  the  acoustic  normal  modes  the  equation  to  be  solved  is  the 
Helmholtz  Equation 

d2Zn        2 


+4,^  =  0  (i) 


2       '  '^zn^n 


where 

2 
2         _CO 2 

Kzn  ~~       j        K" 
C 

is  the  vertical  wavenumber,  Z„  is  the  n'h  normal  mode,  z  is  the  vertical  coordinate,  co  the 
acoustic  angular  frequency,  c=  c(z)  the  sound  speed  and  k„  the  horizontal  wavenumber. 
With  a  pressure  release  surface  and  a  rigid  bottom  the  appropriate  boundary  conditions 
are 

Zn  =  0  at  the  surface  (2) 

dZn 

— ; —  =  0  at  the  bottom    .  (3) 

dz 

Normal  modes  subjected  to  general  boundary  conditions  are  discussed  in  Appendix  B. 
Note  that,  in  (1)  through  (3),  we  have  supressed  the  dependence  in  range  which  has  no 
effect  on  the  local  modes.    Equation  (1),  together  with  (2)  and  (3),  is  a  Sturm- Liouville 


problem  where  the  Zn's  are  the  eigenfunctions  and  the  k„'s  are  the  eigenvalues.    The 
eigenfunctions  can  be  normalized  using  a  normalization  constant  C„  such  that 

Zn  =  CnZn  (4) 

where 

jZnZmdz  =  6nm    .  (5) 

Note  that  the  integration  is  over  the  entire  depth  in  (5),  5nm  is  the  Kronecker  Delta,  and 
Zn  and  Zm  are  orthogonal  for  n  #  m. 

There  are  several  ways  to  obtain  the  first  order  WKB  solution.  A  physical  ap- 
proach is  to  consider  the  transmission  of  plane  waves  through  a  layered  medium.  The 
WKB  solution  is  obtained  as  we  let  the  the  layer  thickness  approach  zero.  This  physical 
approach  will  be  discussed  next. 
2.     Physical  Approach 

In  each  isospeed  layer  the  solution  to  the  Helmholtz  Equation  is  given  by 

Zn  oc  e     m 

where  K'2n  is  the  vertical  wavenumber  in  the  j'h  layer.    The  transmission  coefficient  from 
layer  j  to  layer  j  +  1  is  given  by  (Kinsler  et  al.,  19S2) 

A.  -»« 


or 


Kzn~i~Kzn 


r„.,= 1 


>J+'  A^„ 
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2     kL 


Akj 
where  Ak{„  =  k{;x-k^.   With  — j—  <^1  we  get 


1    AkJb 


By  letting  |  Z\  \  =  1  we  have 


/-I 

n  rm,  m+,  w>* .  «-i  x^^ 


where  Z^  is  the  value  at  the  interface  of  the/*  layer  and  the  (/'  +  l)'*  layer  and  Azm  is  the 
thickness  of  the  mth  layer.   Taking  the  thickness  of  the  layers  to  zero  we  have 


1    ['mdKm        f« 


or 


where  kc  is  the  vertical  wavenumber  at  z  =  z,  and 


-J''- 


0  =  |  K,„  |  <fe      . 


The  vertical  wavenumber  Km  generally  can  take  on  real  or  imaginary  values.    If 
k„  is  real,  i.e.,  k^  >  0,  the  solution  is 

'    i<t>  i  u>    —'<P 
„  .  ,       a  e    +b e 

Zn(z)  = 7= 

\JKzn 

or 

a  cos  0+6  sine/) 
A,(z)  = = (°) 

V  Kzn 

where  a'  and  b'  or  a  and  £  are  constants  to  be  determined  by  normalization  and  the 
boundary  conditions.   On  the  other  hand  if  k„  is  purely  imaginary,  i.e.,  k\„  <  0,  we  have 

n^)  =  ~~fT=T~     ' 

V  I  Kzn  I 


These  last  two  formulae  are  precisely  the  WKB  solution  to  the  Helmholtz  Equation  in 
regions  not  close  to  a  point  where  k2„  =  0  (i.e.,  a  turning  point). 
3.     Basic  Formulae 

By  substituting  the  WKB  solutions  (6)  or  (7)  in  the  Helmholtz  Equation,  one 
can  easily  show  that  they  are  exact  solutions  if 

Therefore,  if  |r(z)|<£l  the  WKB  solutions  (6)  and  (7)  are  good  approximations  to  the 
exact  solution  for  k\„  >  0  and  k£,  <  0,  respectively.  In  Appendix  A  a  more  detailed  and 
mathematical  derivation  of  the  WKB  solution  and  its  validity  is  presented.  Let  us  call 
these  solutions  the  first  order  WKB  approximation.  But  (6)  and  (7)  are  not  valid  when 
k\„  =  0  (i.  e.,  at  a  turning  point)  or  even  in  a  region  where  k\„  —  0.  We  need  a  solution 
valid  at  and  near  the  turning  point  (i.e.,  in  the  critical  region)  to  make  the  liaison  be- 
tween the  oscillatory  region  (k\^>0)  and  the  exponential  region  (k^O). 

Let  us  suppose  that  there  is  a  turning  point  at  z  =  0  and  k\„  >  0  for  z  >  0  and 
kI,  <  0  for  z  <  0.  Consistent  with  a  first  order  approximation,  let  us  assume  that  near  the 
turning  point  k]„  =  yz  where 


With  a  change  of  coordinate 


the  Helmholtz  Equation  becomes 


C  =  -y1/3z   , 


d2Zn 

--CZ  =0 


which  is  the  Airy  Equation  with  the  general  solution  given  by 

Z„(0  =  aAi(0+bBi(Q    . 
The  Airy  Functions  can  be  expressed  as,  with  t,  >  0  , 


Ai(-Q  = 


v  - 


2    ,3/2 


-i/3iy<. 


+J< 


2    r3/2 


/3\     3 


^-0-7fp-/3(f^2)-y"<T^'2) 


/f«)  = 


Bi(Q  = 


,'C 


''-„(fc-)-/,;3(j^ 


I?3'2  Wife"2 


where  7's  are  Bessel  Functions  and  /'s  are  Modified  Bessel  Functions. 

Now  we  have  a  complet  set  of  first  order  approximate  solutions  covering  the 
entire  range  of  x\n.   Summarizing,  we  have: 

(a)  in  the  oscillatory  region  (k^>0) 


Zn{2)  - 


a,  sin  </>+&i  cos  0 


(8) 


(b)  in  the  critical  region  (kI„  ~  0) 


Kln  =  vz 


1/3 

C  =  -y      Z 


Zn(z)  =  a2Ai(Q+b2BiC)    , 


(9) 


(c)  in  the  exponential  region  (k*„<0) 


Zn{z)  = 


a3e  *+b^ 


V  I  Kzn  I 


(10) 


To  assure  continuity  in  Zn{z)  at  the  boundaries  between  regions,  we  must  relate 
the  constants  a,  and  b,  carefully.  An  asymptotic  expansion  of  the  critical  region  solution 
for  large  positive  values  of  K\n  (or  — C  >  1)  must  match  the  solution  in  the  oscillatory  re- 
gion and  for  large  negative  values  of  k]„  (or  C  >  1)  must  match  the  solution  in  the  expo- 
nential region.   The  connection  formulae  between  the  three  regions  are  given  by 


V  i  Kzn 


oAi{Q    -»      sin(0+-Y) 

\  Kzn 


oBi{Q 


^JK2 


cos(</)+  — ) 


where  the  left  sides  match  the  Airy  Functions  for  k^O  and  the  right  sides  match  for 
k*>0  (Bender  and  Orszag,  1978)  with 


o  =  yiny 


-1/6 


Therefore,  (8),  (9)  and  (10)  can  be  recast  respectively  as 
(a)  oscillatory  region  solution 


Z„(z)  = 


a  sin(0+  —  )+b  cos(0+  — ) 
m  -4 


(ll.a) 


(b)  critical  region  solution 

Zn(z)  =  acjAi(0+boBi(0    , 


(c)  exponential  region  solution 


Z„(z)  = 


(al2)e~<p+be' 


j\Kzn 


(11.6) 


(He) 


The  constants  k„,  a  and  &  are  determined  by  normalization  and  the  boundary 
conditions.  The  equations  for  the  eigenvalues  k„,  i.e.,  the  characteristic  equations,  in- 
volve the  solution  in  the  oscillatory  region  (ll.a),  as  will  be  discussed  later. 

Equations  (ll.a,b,c)  are  not  easy  to  use  in  the  computation  of  normal  mode 
shapes.  Where  does  one  stop  with  one  formula  and  start  with  another?  This  problem  is 
avoided  if,  after  the  determination  of  the  eigenvalue  k„  and  the  constants  a  and  b,  the 
computation  of  the  normal  mode  is  done  using  the  Langer  Formula  (Nayfeh,  1973; 
Bender  and  Orszag,  1978) 


Zn(z)  =  2N/*(y  101 


1/6 


1 


aAi 


i  -4V 


+bBi 


<h) 


2/3 


(12) 


with 


</)  =       |  Kzn  \dz 


2n 

0 


and 

s  =  -zl\z\    . 

This  formula  has  the  advantage  of  giving  a  single  and  continuous  solution  for  the  entire 
range  of  k\„.  Bender  and  Orszag  (1978)  have  shown  that  for  k\„  ~  0  and  K2m  =  yz  the 
formula  gives  exactly  (1  l.b),  for  k£>0  it  asymptotically  approaches  (11. a)  and  for  k\n<^Q 
it  approaches  (1  l.c). 

In  rectrospect,  in  order  to  obtain  the  first  order  WKB  solution,  an  algorithm 
must  include  the  following  steps: 

(a)  application  of  a  boundary  condition  to  (ll.a,b,c)  to  get  the  relationship 
between  the  constants  a  and  b; 

(b)  application  of  the  other  boundary  condition  to  (11. a)  to  get  the  character- 
istic equation; 

(c)  with  (12)  compute  Z„{z); 

(d)  with  (4)  and  (5)  normalize  Zn(z). 
4.     Comments 

The  WKB  solution  is  a  high  frequency  approximation.  This  means  that  it  gets 
better  as  the  frequency  gets  higher.  It  must  be  noted  that,  although  Normal  Modes  are 
a  full-wave  exact  solution  to  the  wave  equation,  the  WKB  modes  are  approximate  sol- 
utions and  their  accuracy  is  frequency  dependent. 

B.     DETERMINATION  OF  WKB  MODE  PARAMETERS 

Now  that  we  have  a  uniformly-valid  first  order  solution  to  the  Helmholtz  Equation, 
we  are  ready  to  apply  the  boundary  conditions  to  get  expressions  for  the  k„'s  and  the 
constants  a  and  b  (or  their  equivalents).  There  are  four  classes  of  normal  modes  to 
consider.  Each  class  has  different  mathematical  expressions  for  the  mode  parameters  (a, 
b  ,  or  their  equivalents,  and  k„).  Therefore,  for  an  arbitrary  sound  speed  profile  the  first 
procedure  for  normal  mode  calculation  is  to  determine  which  class  each  mode  falls  into 
and  then  go  through  the  steps  described  at  the  end  of  the  previous  section. 
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1.     General  WKB  Formulae 

We  start  by  deriving  the  more  general  formulae  which  are  applicable  to  modes 
that  do  not  have  turning  points  in   close  vicinity  of  the  boundaries. 
a.     Class  I:  Pressure  Release  -  Turning  Point 

First  let  us  consider  a  mode  whose  oscillatory  region  is  bounded  by  a  turn- 
ing point  at  the  depth  z  —  z  and  the  surface  (z  =  0).  With  the  bottom  at  z  =  //,  the 
boundary  conditions  are 

Z„(0)  =  0  (13) 

dZn 
~^(H)  =  0    .  (14) 

We  will  call  this  class  of  modes  as  PR-IP  (Pressure  Release  -  Turning  Point). 

In  the  exponential  region  (z4z  <  H),  the  WKB  solution  can  be  recast  as 

cosh(  4>-4>b) 
Z„oc (15) 

V  I  Kzn  I 

with 

<f>  =  I    |  K2n  I  dz. 

Jz 

In  order  to  satisfy  the  boundary  condition  (14),  4>b  must  be  given  by 

*»-!" k'"Uz~tanh"'(2kli2  rfl£"' )—  ' 

We  need  to  connect  (15)  to  the  solutions  in  the  critical  and  oscillatory  regions  next.  After 
application  of  the  connection  formulae,  the  results  in  the  three  regions  are: 

(a)  exponential  region 

Zn  =     i =-  cosh(<t>-<j>b)    ,  (16.a) 

V  I  K2  n  I 

(b)  critical  region 
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C  =  -y1/3(z-z) 


Zn  =  e%^/(0+-~-^'(0    , 


(16.*) 


(c)  oscillatory  region 


Z*~ 


sin(0+  -|- )+     gJl  cos(<£+  -J- ) 


2Vk 


(16.c) 


with 


<t>  =\    \K2n\dz 


in   this   region.      The   corresponding   Langer   Formula  that   asymptotically  matches 
(16.a,b,c)  is 


Z„  =  2V,T^U|V'6-J=  lt*Ai 


— V-      3 


VK 


|^ 


[|0V/3 


>(17) 


with 


5  = 


A 

z—z 


z-z 


and 


<P 


-JS 


rfz. 


As  the  oscillatory  region  solution  must  satisfy  the  surface  boundary  condi- 
tion, we  have 


.-** 


2A 


cos 


)>. 


A 


sin 


):• 


™^z+  "J1 


=  0    , 


(18) 


or 
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1   \  ..-1/    e 


K,„dz  =  [  n — r  \n—  tan 


-4> 


v         4    /  \    ->,*„ 


o  \  /  \  ^ 

This  is  the  characteristic  equation  for  PR  -  TP  normal  modes.  The  solutions  to  the 
characteristic  equation  give  the  eigenvalues  k„'s. 

b.     Class  II:  Turning  Point  -  Rigid  Bottom 

The  next  class  of  normal  modes  to  be  considered  has  the  oscillatory  region 
between  a  turning  point  and  the  bottom.  It  will  be  called  TP-RB  (Turning  Point  -  Rigid 
Bottom).  For  this  class,  let  us  express  the  solution  in  the  exponential  region,  0<z<;i, 
as 

Zn  oc  — ,  sinh((/>— 4>s) 

\l  I  Km  I 

with 


<t>  —  J      I  Km  I  dl 


and 


?2 


4>s  = 


K.n\dz. 


zn 

0 


Using  the  connection  formulae,  we  get  for  the  three  regions: 
(a)  Exponential  Region 

Zn  = r. sinh(<b-<bs) 


(b)  Critical  Region 


V  I  Km  I 


2  /        A\ 


C  =  -y1/V^A) 


Zn  =  e*<oAm—2-^-aBm 
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(c)  Oscillatory  Region 


Zn  = 


4rsin(^+-S-)- 


-<P. 


2yJ  Kzn 


COS 


(*+i) 


(19) 


with 


4>=j\k2„\ 


dz 


in  this  region  (z<£z  <  H).   The  appropriate  Langer  Formula  for  this  class  is 


Z„  =  2V*(  y  |0| 


2.\4,\\l*-J=\e*.Ai 


[|0V/3 


,"*i 


5/ 


(|  0)^)1(20) 


where 


z—z 


z-z 


and 


0 


-J>- 


|</z. 


Equation  (19)  must  satisfy  the  bottom  boundary  condition  (14),  therefore,  we  have 


sin 


T-- 


r«^+y+tan  '| 


<A+  A  e~*.  s,  -, 


-<ps 


■ZV 


=  0 


(21) 


where 


Z)  = 


_1__^£L  \ 

2kL  ^  rH 


(22) 


It  follows  from  (21)  that  the  characteristic  equation  for  TP  -  RB  is 
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K,„dz  =  (  n——r  )n—  tan 


j  K„rfz -(«-t) 


^-Z><< 


for  Z)  <  0,  and 


Ja  >wfe=  U-  — JTr-tan" 


e*<+  -f  <T*< 


>  ±^—De*' 


for  Z)  >  0. 

c.     Class  III:  Turning  Point  •  Turning  Point 

The  third  class  of  modes  has  the  oscillatory  region  between  two  turning 
points,  at  z  =  z,  and  z  =  z2  with  z,  <  z2.  They  will  be  refered  as  TP-TP  (Turning  Point  - 
Turning  Point). 

Let  us  express  the  solution  in  the  exponential  region  near  the  surface, 
0<,z4z{,  as 

c, 
Zn  = ,  sinh(0 ,  -(f) s) 

V  I  Kzn  I 

where 


4>\  =\     \*m\dz 

Jz 


and 


^5=       \f<zn\dz, 
Jo 

and  the  solution  in  the  exponential  region  near  the  bottom,  z2<£z  <  H,  as 


Zn  =     j==  cosh[4)2 -<f>b) 
V  I  Kzn  I 

with 
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02"  =  J^   \Kin\dz 


and 


In  addition,  let  us  express  the  solution  in  the  critical  region  centered  at  z  =  z,  as 

Z„  =  Cle*'«rtXi(Cf)-  -y  <?"%/%)  (23) 


with 


and 


.2      _  /  A    >. 


C,  =  -yj'V^i). 


and  in  the  other  critical  region,  centered  at  z  =  z2,  as 

Z*  =  c2eK2AiC2)+  4f-  e-+>o2Bi&2)  (24) 


with 


and 


4»  =  y2&-z)    . 


f~"   -1/6 
°2  =  \ln  Vl 


r  1/3/ A         \ 

C2  =  -y2'  (z2-z). 


By  letting 
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1 

c,  = — , 


>+^ 


and  connecting  (23)  to  the  solution  in  the  oscillatory  region  zx<^z<^z2,  we  get 

Zw  =  ^i=-sin(^+-|— «i)  (25) 

V  Kzn 

where 


*?-/> 


,i^Z 


and 


<-*• 


"1  =  tan  V"^/ 

Note  that,  in  obtaining  (25),  the  trigonometric  identity 


a:+£2  sin(Wtan  'y)  (26) 


has  been  used.    In  the  same  way,  by  letting 


c->  = 


tf  ^+ 


and  connecting  (24)  to  the  oscillatory  region  solution  we  get 


Zn  =  -4=-sln(fc+f+a2)  (27) 

\lKzn 


where 


02  =  I     \xzn\dz 

"2 
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and 

a,  =  tan       —    . 

V  2e^  J 

Realizing  that  the  two  solutions,  (25)  and  (27),  in  the  oscillatory  region 
must  be  identical,  we  obtain 

sin(0t+^--«i)  =  sin(07+^-+a2)    .  (28) 

Using  (28),  the  characteristic  equation  for  this  class  of  normal  modes  is  obtained: 


Ja  K2ndz=  ffl-yj)i+ar«; 


To  compute  the  normal  mode  we  can  use  (20)  for  z  <,  z2  and  (17)  for  z  >  z2. 
d.     Class  IV:  Pressure  Release  -  Rigid  Bottom 

Finally,  we  must  consider  modes  having  no  turning  points.  We  call  this  class 
PR-RB  (Pressure  Release  -  Rigid  Bottom).  Here  we  express  the  solution,  which  is  always 
oscillator)',  as 

Zn  oc —j=- sm  <f)-\ — ■ cos  <f)  (29) 

V  Kzn  V  Kzn 


with 


J0 


<t>=\    \f<2n\dz    . 

Since  the  surface  boundary  condition  must  be  satisfied,  we  must  have  b  =  0,  and  hence 

Z„  =  -i=- sin  </>    .  (30) 

On  the  other  hand,  (30)  must  satisfy  the  bottom  boundary  condition  (14)  and  this  leads 
to  the  characteristic  equation  equation  for  the  PR  -  RB  modes: 
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?H 


K2ndz  =  nn—  tan     (  — 


for  D  >  0,  and 


i</z  =  nn+  tan 


-l 


i) 


for  D  <  0,  where  D  is  given  in  (22). 

2.     Formulae  Associated  With  Near-Boundary  Turning  Points 

In  our  derivation  of  the  formulae  in  the  previous  section,  we  have  applied  the 
boundary  conditions  to  the  oscillatory  or  exponential  region  solutions,  assuming  that  the 
boundaries  are  nowhere  close  to  any  turning  point.  In  the  special  case  that  a  boundary 
lies  inside  a  critical  region,  the  respective  boundary  condition  must  be  applied  to  the 
critical  region  solution  instead.  Different  formulae  associated  with  the  first  three  classes 
of  modes  for  this  special  case  will  be  derived  next. 

a.     Class  I 

PR  -  TP  modes  have  lower  turning  points.   Since  the  solution  in  the  critical 
region  must  now  satisfy  the  bottom  boundary  condition  we  obtain 


where 


Zn  =  Bi'ttH)oAi(Q-Ai'ttH)oBi{Q 


iH  =  yXI\H-z) 


and 


y  =  - 


dz 


Ai'  and  Bi'  are  the  derivatives  of  the  Airy  Functions  with  respect  to  £  and  are  defined 
by  (with  C  >  0) 


Ai\-Q  =  -\c 


J 


-2/3 


i^K(-k'2 
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Bi'(-Q  =  -L-l 
-s/3 

^/'(0  =  -yC 

5/'(C)  =  -=-C 
V3     L 

/ 

-2/n  3 


/_,„(  ^c3/2i 


2/3 


i^y 


L2;3(|c"2)-/2/3(|c 


3/2 


-2/3(  3 


C3/2)+/2/3(  -rC 


2/n  3 


_2_v3/2 


and 


Defining 


^i  =  Bi'{U 


Bx  =  Ai'(ZH) 


and  applying  the  connection  formulae,  the  solution  in  the  oscillatory  region  becomes 


.4 


/?, 


Zn  =     _L_  sin(c/>+  -j  ) p==-  cos((/>+  -j  ) 


V  ^Zrt 


An  application  of  the  surface  boundary  condition  gives  the  following  characteristic 
equation: 


\\ndz  =  ("-|)*+  tan"1  (A 


A.     C/ass  // 

TP  -  RB  modes  have  upper  turning  points.  To  satisfy  the  surface  boundary 
condition,  we  require  the  solution  in  the  critical  region  to  vanish  at  the  surface,  i.  e., 


with 


Zn  -  Bi^s)oAi(Q-Ai^s)oBi{Q 


y  1/3  a 


and 


20 


dKi 


y- 


Defining 


and 


A2  -  BiUs) 


B2  =  Ai{Cs) 

and  applying  the  connection  formulae,  the  solution  in  the  oscillatory  region  becomes 

A,       .    /  .      it  \         B, 


Z„  = 


^JKzn 


sin(0+  -S-  )-  -p=-  cos(</>  +  -j  ) 


V  ^;n 


The  subsequent  application  of  the  bottom  boundary  condition  leads  to  the  following 
characteristic  equation  : 


H       .       (        1  \       t     -x(  A2+DB2 
KK,Jz-\n-T)x-X*n    [-^ZoT2 


for  D  <  0,  and 


// 


K7„dz  =  [  n—  —  \n—  tan 


A2+DB2 
B2-DA2 


for  D>0,  where  D  is  defined  in  (22). 
c.     Class  III 

With  the  above  results,  we  can  easily  that  only  two  modifications  in  the 
general  formulae  for  TP  -  TP  modes  are  required.  These  include  replacing  e*>  by  A2  and 
—r—  by  B2,  if  the  upper  turning  point  is  close  to  the  surface,  and  <?**  by  A{  and  — 
by  Bu  if  the  lower  turning  point  is  near  the  bottom. 
3.     Criterion  For  Formulae  Selection 

We  must  define  a  criterion  for  switching  from  the  general  formulae  to  the  for- 
mulae associated  with  near-boundary  turning  points.  It  was  found  by  Bender  and 
Orszag  (1978)  that  the  critical  region  extends  on  each  side  of  the  turning  point  until 
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ICI  ~  1  .In  accordance,  we  will  use  the  general  formulae  unless  the  distance  between  a 
turning  point  and  a  boundary  corresponds  to  values  of  |£|  smaller  than  1. 
4.     Normalization 

Once  normal  mode  parameters  are  computed  using  the  appropriate  formulae 
including  the  characteristic  equations,  the  normalization  of  the  normal  modes  can  be 
achieved  by  numerical  integration  over  depth.  The  normalized  modes  (Z„'s)  are  related 
to  the  Z„  's  by 


z  =c  z 


with 


f*H 


1 
C2 


Z2ndz 


5.     Parameter  Storage 

We  now  define  the  necessary  parameters  required  to  completely  characterize 
each  mode.  The  first  parameter  is  obviously  the  class  number.  By  checking  the  formulae 
developed  in  the  previous  sections  it  is  seen  that  all  the  constants  (0„  4>b,  D,  etc.)  can 
be  computed  from  the  horizontal  wavenumber,  the  depths  of  the  turning  points  and  the 
depth  of  the  ocean.  Therefore,  the  parameters  required  to  parameterize  each  class  of 
modes  are: 

a.  Classes  I  And  II 

1.  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 

4.  Depth  of  the  turning  point 

b.  Class  III 

1.  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 

4.  Depth  of  the  upper  turning  point 

5.  Depth  of  the  lower  turning  point 
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c.     Class  IV 

1.  Class  number 

2.  Horizontal  wavenumber 

3.  Normalization  constant 
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III.     ACCURACY  TEST 

A.     EXACT  ANALYTICAL  SOLUTIONS 

To  quantify  the  accuracy  of  the  first  order  WKB  approximation,  we  will  compare 
the  WKB  solutions  to  exact  analytical  solutions  to  the  Helmholtz  Equation  for  three 
abstract  sound  speed  profiles.  The  exact  solutions  are  derived  in  this  section.    We  use 
as  boundary  conditions  pressure  release  surface  and  rigid  bottom  for  all  three  cases. 
1.     Positive  Exponential  Profile 

In  this  upward  refractive  profile,  sound  speed  increases  exponentially  with  the 
depth  and  is  given  by 

c(z)  =  cQePz 

where  /?  is  a  constant.   The  surface  is  at  z  =  0  and  the  bottom  at  z  =  H.  The  Helmholtz 
Equation  is 


d2Zn 
dz2 


+ 


O) 


Lv  co 


cne 


§2 


2        2 
-Kn 


z„  =  o 


or 


d2Z, 


dz' 


+(Ky^-Ki)zn=o 


(31) 


with  kI  = 


(O2 


°~      r2 


With  a  change  of  variable 


X  =  e 


(31)  becomes 


^i(^)+(^-^»=° 


After  division  by  (]2x2  and  using 


(32) 
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and 


(3?) 


we  obtain 

1     d   (_   dZn  \  Y.2     f? 


x  ^  r  dx    H^-—  )z«  =  0 


This  is  the  Bessel  Equation  and  its  general  solution  is 

Zn  =  aJan(«Qx)+bYan(<x0x)    .  (34) 

The  boundary  conditions  are 

Zn(z  =  0)  =  Zn(X=\)  =  0  (35) 

dZ„  dZ„ 

-^(2  =  H)  =  -^-(X  =  Xh)  =  0    ■  (36) 

The  application  of  (35)  and  (36)  to  (34)  results  in  two  algebraic  equations 

aJ^+bY^-O  (37) 

^'4aoXH)+bY'^0xH)  =  0    .  (38) 

As  this  linear  system  of  algebraic  equations  is  homogeneous,  there  is  a  nontrivial  sol- 
ution only  if  the  corresponding  Jacobian  is  zero,  i.e., 

•V«o)  Y'4«aH)-  Y^oV4aoXn)  =  0    • 

This  characteristic  equation  must  be  solved  in  order  to  obtain  the  a„'s,  which  are  the 
scaled  eigenvalues.  After  the  a„'s  are  evaluated,  the  k„'s  can  be  determined  using  (33).  It 
follows  from  the  surface  boundary  condition,  expressed  in  (37),  that  the  specific  solution, 
aside  from  a  multiplicative  constant,  is 

Zn(z)  oc  rjaoVjao^-yjo^yJotoe-^    . 

The  constant  of  proportionality  is  given  by  normalization. 
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2.     Negative  Exponential  Profile 

The  second  profile  considered  is  a  downward  refractive  one.   Here,  sound  speed 
decreases  exponentially  with  depth  and  is  given  by 

c{z)  =  c0e  F     . 
Thus,  the  Helmholtz  Equation  is 

^+(Kyz-Kl)Zn  =  0  (39) 


d2Z„ 


dz2 

where  k0  =  — . 
With 

and  a0  and  a„  given  by  (32)  and  (33)  respectively,  (39)  can  be  recast  as 

1  "  ^4-4V-o  ■  W 


j  <m  y  <n  j  \ "   ? 

The  general  solution  to  (40)  is 

Zn  =  aJ,(a0£)+bY3n(v.0Z)    . 

Using  the  same  anology  as  in  the  previous  case,  we  obtain  as  characteristic 
equation 

^WW^^vMb)  =  o  (41) 

and  specific  solution 

z„(2)ocr3n(ao)^(«0/7)-^n(oto)^n(«o^)  • 

The  k„'s  can  be  found  using  (33)  after  solving  (41)  for  the  a„'s. 
3.     Hyperbolic  Profile 

In  this  third  case  we  assume  that  the  sound  speed  has  a  minimum  at  z  =  z0.  The 
sound  speed  profile  is 

c(z')  =  c0  cosh(/?z') 
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with 


z   =  z— ; 


The  Helmholtz  Equation  to  be  solved  is  now 


SzK 


dz 


,2 


+ 


(jo_y — ! 

v  c°  >    cosh' 


W) 


—K, 


z„  =  o 


(42) 


With  the  following  change  of  coordinates 

X  =  tanh(/fe') 


(42)  becomes 

K-fi-k 

■>    dZn 

+[K20(l-X2)-Kl]Zn  =  0 


(43) 


After  dividing  (43)  by  P2{l~x2)  and  using  the  definitions  given  in  (32)  and  (33),  (43)  can 
be  recast  as 


d 
dx 


{X~X)~ 


+H-7^)Z"  =  0 


(44) 


This  is  the  Associated  Legendre  Equation.  The  general  solution  to  (44)  is 

Zn  =  <(x)+bQ:(x)    , 
with 

u  =  a„ 


and 


v(v+l)  =  a0 


The  functions  P»Xx)  an<^  QXx)  are  the  Associated  Legendre  Functions  of  the  first 
and  second  kind  of  order  fj.  and  degree  v  .  Applying  the  boundary  conditions  we  get  the 
system  of  equations 


<\xs)+t>Q:(xs)  =  o 


(45) 
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with 

XS  =  tanh(-/?z0) 

and 

Z//  =  tanh[/?(tf-z0)]    . 
Thus,  the  corresponding  characteristic  equation  is 

p:(xs)q':(xh)-q:(xs)p":(xh)=q  ■  w 

Solving  (46)  for  n  gives  the  eigenvalues  a„  's  and  hence  the  horizontal  wavenumbers 
■k„  's.   Using  (45)  the  solution  becomes 

Zn{z)  oc  Ql»{xs)Kn(  tanh  Pz')-Pav«{Xs)Q*A  tanh  ffl    . 

B.     ANALYTICAL  EVALUATION  OF  WKB  PHASE  INTEGRALS  AND 
DERIVATIVES 

Now  we  use  the  WKB  formulae  to  obtain  the  first  order  WKB  solution  for  the  three 
analytical  sound  speed  profiles.  The  objective  in  this  section  is  to  evaluate  analytically 
all  the  related  WKB  integrals  and  derivatives  so  that  numerical  errors  can  be  eliminated 
in  one  of  the  error  analysis.  Horizontal  wavenumbers  computed  using  both  the  exact  and 
approximate  (i.e.,  WKB)  characteristic  equations  will  be  compared  in  the  next  section. 
The  evaluation  of  the  horizontal  wavenumbers  is  emphasized  because  they  are  the  key 
parameters  in  normal  mode  computations. 
1.     Positive  Exponential  Profile 

The  vertical  coordinate  z  is  positive  downward,  the  surface  is  at  z  =  0  ,  the  bot- 
tom at  z  =  II  and  the  sound  speed  profile  is 

c{z)  =  cQeH     . 

With 

the  corresponding  wavenumber  profile  is 
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Let  us  first  consider  modes  that  have  a  turning  point  at 
(oscillatory  region),  kz„  is  given  by 


z.    TorO<z<z 


/    -202       2 

K2n  ~  Ko\ e       y» 


with 


Vn  Kr 


The  spatial  phase  4>  in  this  region  is  given  by 


</>  = 


1*2 


I  K2n  |  dz 


(47) 


With  the  change  of  coordinate 


(47)  becomes 
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In  the  region  z  <  z  <,  II  (exponential  region),  we  must  use 


/    2        2-202 


I  »cr„  I  =  v'  Kn-K0e 


The  spatial  phase  is  now 


or,  with 


4>  -        I  **«  I  <&     , 
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We  must  also  compute 


D  = 


c!\k. 


2\k2J 


dz 


7=tf 


which  can  be  equated  as 


D  = 


4fi 


,-fiH 


2      (4-4e-V")312 


Figure  1  shows  a  PR  -  TP  mode  in  this  upward  refractive  medium. 
If  the  normal  mode  is  PR-RB  the  phase  is  given  by 
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dz 


or 


4>  = 
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and 


Z)  =  - 
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For  each  case  we  must  also  apply  the  respective  characteristic  equation. 
2.     Negative  Exponential  Profile 

For  this  profile  we  have  two  types  of  normal  modes:  TP  -  RB  and  PR  -  RB. 
The  sound  speed  is 
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Figure   1.       PR  -  TP  Mode 


c{z)  =  c0e  p2 


Consequently,  for  a  TP  -  RB  normal  mode  and  for  z  >  z  (oscillatory  region)  we  have 


\KzJ  =KQ\leiP2~y1n 


So  in  this  region  the  phase  is 
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4> 
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with 
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Also  we  have 


Z)  = 
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2/?H 


2      (k^"-k^2 


(48) 


In  the  exponential  region  (0  <  z  <  z), 


'  i       ipz 
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and  thus  the  phase  is 
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Figure  2  shows  a  TP  -  RB  mode  in  such  a  downward  refractive  medium. 
For  a  PR  -  RB  normal  mode  we  have 


4>  = 


x-y2-y*tan  '| 


x-vl 


^2 


-1*00 


"W) 


and  D  is  given  by  (48).    Again,  for  each  of  the  classes  the  respective  characteristic 
equation  must  be  used. 
3.     Hyperbolic  Profile 

We  consider  two  types  of  normal  modes:  TP  -  TP  and  PR  -  RB.    For  the  type 
TP  -  TP  and  in  the  oscillator)'  region  (z\<z  <  z'2)  the  phase  will  be 
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Figure  2.      TP  -  RB  Mode 


4>  =\     I  O  dz' 


(49) 


Using 


£  =  tanh(/?z')    , 
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an  = 
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and 


(49)  becomes 
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For  the  exponential  region  near  the  bottom  (H>  z'  >  z'2)  the  phase  is 
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For  this  type  of  normal  mode  D  is  given  by 


D  = 


In  the  exponential  region  near  the  surface  (£',  >  z')  the  phase  is 
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If  the  mode  is  PR  -  RB  the  phase  will  be 
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and  D  is  given  by 
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Figure  3  shows  a  TP  -  TP  mode  trapped  by  the  refractive  index. 

C.     COMPARISON  OF  RESULTS 

1.     Horizontal  Wavenumber  Error  Analysis 

The  exact  characteristic  equations  derived  in  section  A  for  the  three  abstract 
profiles  were  solved  using  iterative  procedures  to  obtain  the  benchmark  Kn  's.  Similarly, 
the  approximate  WKB  characteristic  equations  developed  in  the  previous  chapter  were 
also  solved  for  the  three  profiles  with  the  use  of  the  algebraic  equations  developed  above. 
Very  low  frequencies  were  chosen  for  the  comparison  intentionally.  This  would  give  the 
WKB  method  a  real  test,  since  WKB  is  inherently  a  high-frequency  approximation.  For 
the  exponential  profiles,  we  used  co  =  10  s-1  or  a  frequency  of  1.59  Hz.  For  the 
hyperbolic  profile  the  frequency  used  was  even  lower,  it  was  0.23  Hz.  Some  of  the 
computed  horizontal  wavenumbers  (k„)  are  presented  in  Tables  1  through  5.  The  unit  for 
k„  is  inverse  meter  (m_1).  The  depth  of  the  ocean  is  taken  to  be  10000  m. 

In  Tables  1  and  2  the  exact  and  WKB  results,  as  well  as  the  absolute  and  relative 
errors,  for  the  positive  exponential  profile  are  displayed.  Tables  3  and  4  show  the  results 
for  the  negative  exponential  profile.  In  Table  5,  the  exact  and  approximate  results  for 
the  hyperbolic  profile  are  compared.  At  0.23  Hz,  mode  2  is  TP  -  TP  and  mode  3  is  PR 
-  RB. 

Overall,  we  can  see  that  the  absolute  error  (the  absolute  value  of  the  WKB  re- 
sult minus  the  exact  result)  varies  between  10"4  and  10~7  m~x  and  the  relative  error  (the 
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Figure  3.      TP  -  TP  Mode 

absolute  error  divided  by  the  exact  value)  varies  between  10~2  and  10~5.  It  is  noted  that 
for  the  modes  with  turning  points  not  in  close  vicinity  to  the  boundaries  the  absolute 
error  usually  stays  below  10~5  nrl.  The  fact  that  the  error  increases  when  a  turning 
point  is  close  to  a  boundary  is  easily  explained.  The  solutions  in  the  critical  regions  are 
the  less  accurate  because  they  use  a  linear  approximation  for  k\„.  So  when  they  are  used 
to  match  the  boundary  conditions  the  error  increases. 
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Table 

1.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 

Mode 

# 

Exact  k„  {m'1) 

WKB  Kn  (m-1) 

Absolute 
Error  (m_1) 

Relative 
Error 

5 

.43S61  x  10-2 

.43866  x  10-2 

5  x  10-7 

1.1  x  10-4 

6 

.40441  x  10-2 

.40445  x  10  2 

4  x  10-7 

1.0  x  10-4 

7 

.37227  x  10-2 

.37230  x  10-2 

3  x  10-7 

.8  x  10-4 

8 

.34179  x  10-2 

34182  x  10-2 

3  x  10-7 

.9  x  10-4 

9 

.31274  x  10-2 

.31277  x  10-2 

3  x  10-7 

1  x  10-4 

10 

.28556  x  10-2 

.28568  x  10-2 

1.2  x  10-6 

4.2  x  10-4 

Table 

2.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 

Mode 

Exact  Kn  (w_I) 

WKB  k„  {m-') 

Absolute 
Error  (m_I) 

Relative 
Error 

12 

.23190  x  10-2 

.23041  x  10-2 

1.5  x  10  5 

6.3  x  10-3 

13 

.18524  x  10-2 

.18490x  10-2 

3.4  x  10* 

1.8  x  10-3 

14 

.10756  x  10-2 

.10736  x  10-2 

2  x  10~6 

1.9  x  10^3 

Table  3.     NEGATIVE  EXPONENTIAL  PROFILE,  TP  -  RB  MODES 


Mode 

Exact  k„  (m_l) 

WKB  Kn  (m-1) 

Absolute 
Error  (w_1) 

Relative 
Error 

23 

.82307  x  10-2 

.82308  x  10-2 

1  x  10-7 

1  x  10-5 

24 

.79463  x  10-2 

.79464  x  1 0-2 

1  x  10-7 

1  x  10-' 

25 

.76664  x  10-2 

.76665  x  10-2 

1  x  10-7 

1  x  10-5 

26 

.73903  x  10-2 

.73904  x  10-2 

1  x  10-7 

1  x  10-5 

27 

.71 158  x  10-2 

.71155  x  10-2 

3  x  10-7 

4  x  10-5 

In  general,  for  each  type  of  mode  the  WKB  approximation  gets  better  as  mode 
number  increases.  The  only  exception  is  when  turning  points  are  close  to  the  boundaries. 
To  see  this,  let  us  examine  the  results  for  the  positive  exponential  profile  (Tables  1  and 
2).  Starting  from  the  lowest  modes  that  have  one  turning  point,  the  error  decreases  as 
mode  number  increases.  At  mode  number  10,  the  accuracy  of  the  WKB  method  de- 
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Table  4.     NEGATIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

rr 

Exact  k„  (wr1) 

WKB  Kn  (m-1) 

Absolute 
Error  (m_1) 

Relative 
Error 

29 

.65250  x  10~2 

.65091  x  10-2 

1.6  x  10-5 

2.4  x  103 

30 

.61746  x  10-2 

.61683  x  10-2 

1.6  x  10-5 

1.0  x  10-3 

31 

.57723  x  10-2 

.57688  x  10-2 

3.5  x  10-6 

6.1  x  10-4 

32 

.53082  x  10-2 

.53058  x  10-2 

2.4  x  10-6 

4.5  x  10-4 

33 

.47670  x  10-2 

.47652  x  10-2 

1.8  x  10-6 

3.8  x  10^ 

Table  5.     HYPERBOLIC  PROFILE 


Mode 

Exact  k„  (m_1) 

WKB  Kn  (m-1) 

Absolute 
Error  (m_1) 

Relative 
Error 

2 

.901 14  x  10-3 

.S9881  x  103 

2.3  x  10-6 

2.2  x  10-3 

3 

.74139  x  10-3 

.72953  x  10-3 

1.2  x  10-5 

1.6  x  10-2 

creases  because  the  turning  point  comes  too  close  to  the  bottom  boundary.  For  mode 
12  and  higher  a  turning  point  does  not  exist  and  the  oscillatory  region  is  bounded  by  the 
physical  boundaries.  After  the  change  of  mode  class  the  error  starts  to  decrease  again 
as  mode  number  increases. 

Errors  in  the  xn's  for  all  cases  are  small  enough  that  they  do  not  significantly 
influence  the  vertically  integrated  phases  0  because  the  maximum  ocean  depth  is  of  the 
order  of  km's.  In  the  horizontal  direction  the  horizontal  phase  is  approximately  K„r  (this 
phase  is  exact  when  c=  c(z)).  The  error  in  k„  limits  the  range  for  which  the  use  of  WKB 
is  accurate.  If  we  want  to  keep  the  phase  error  below  a  few  degrees,  the  tolerable  error 
in  k„  is  of  the  order  of  10-5  nrx  for  a  range  of  10  km,  10-6  m~l  for  100  km  and  10-7  m'1 
for  1000  km.  All  the  k„'s  associated  with  the  modes  that  do  not  have  interactions  with 
the  bottom  boundary,  as  calculated  from  the  WKB  method,  have  errors  less  than 
10"6  nr\  implying  that  the  method  is  good  to  at  least  100  km  at  1  Hz.  For  higher  fre- 
quencies the  results  would  be  much  better. 
2.     Errors  In  The  Interference  Distances 

The  acoustic  pressure  square  at  far  field  can  be  written  as  (assuming  no  range 
dependence) 
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oo  2£^    ^ 
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where  the  An's  are  constants  and  the  acoustic  source  is  at  z  =  z0  and  r  =  0  (Clay  and 
Medwin,  1977).  We  can  see  that  the  interference  terms  are  functions  of  the  differences 
of  horizontal  wavenumbers.  The  interference  distance  (or  interference  wavelength)  de- 
fined bv 


''/im 


In 


is  therefore  an  important  parameter  for  transmission  loss  calculation  (Chiu  and  Ehret, 
1990).  In  general,  dominant  interferences  are  between  adjacent  modes  (Chiu  and  Ehret, 
1990).  So  let  us  see  what  is  the  size  of  the  errors  in  Ak„  =  \k„+1— k„\  computed  using  the 
WKB  approximation.  In  Tables  6  and  7  we  display  the  exact  and  WKB  Ak„'s,  as  well 
as  the  absolute  errors.  As  expected,  we  can  see  that  the  errors  in  Ak„  van'  in  a  similar 
way  as  in  k„  and  are  slightly  smaller. 

Table  6.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 


n 

Exact  Ak„  {m~l) 

WKB  Ak„  (m-1) 

Absolute  Error 
(m-») 

5 

.3420  x  103 

.3421  x  10-3 

1  x  10-7 

6 

.3214  x  10~3 

.3215  x  10-3 

1  x  W'7 

7 

.3048  x  10-3 

.3048  x  10-3 

~  0 

8 

.2905  x  10-3 

.2905  x  10^3 

-0 

9 

.2718  x  10-3 

.2709  x  10~3 

9  x  10-7 

3.     General  Numerical  Method 

The  WKB  results  in  the  previous  section  were  obtained  using  analytical  means 
rather  than  numerical  methods.  The  reason  for  that  was  we  wanted  to  see  how  accurate 
the  WKB  approximation  is  regardless  of  the  numerical  methods  used  to  evaluate  inte- 
grals and  derivatives.  As  we  see,  absolute  errors  in  the  order  of  10-7  rrrx  for  the  k„'s  are 
achieved  at  1  Hz.  This  means  that  the  WKB  solution  is  very  accurate,  especially  because 
we  used  very  low  frequencies. 
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Table 

7.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 

n 

Exact  Ak„  (w_1) 

WKB  Ak„  (w-1) 

Absolute  Error 
(m-1) 

12 

.4666  x  10-3 

.4551  x  10-3 

1.15  x  10  5 

13 

.7768  x  10-3 

.7754  x  10-3 

1.4  x  10-6 

For  an  arbitrary  profile  all  the  integrations  and  differentiations  must  be  done 
numerically.  It  is  expected  that  the  errors  will  increase  due  to  numerical  noise.  The 
Fortran  program  WKBGEN  (Appendix  C)  developed  in  this  thesis  can  be  applied  to  an 
arbitrary  profile.  For  each  mode  it  finds  class,  horizontal  wavenumber  and  depths  of  the 
turning  points.  This  program  was  also  applied  to  the  three  abstract  profiles.  Some  nu- 
merical results  for  the  k„'s  are  presented  in  Tables  8  through  11. 

Tables  8  and  9  show  the  exact  and  numerical  WKB  results  and  the  respective 
absolute  errors  for  the  positive  exponential  profile.  Tables  10  and  11  show  to  the  cor- 
responding results  for  the  negative  exponential  profile. 

Table  8.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  TP  MODES 


Mode 

u 

Exact  Kn  (/7T1) 

WKB  Kn  (m-») 

Absolute 
Error  (m_l) 

5 

.43S61  x  10-2 

.43874  x  lO"2 

1.3  x  10-6 

6 

.40441  x  10~2 

.40454  x  Ur2 

1.3  x  10-6 

7 

.37227  x  lO"2 

.37234  x  10~2 

7x  10-" 

8 

.34179  x  10'2 

.34184  x  L0-2 

5  x  10-7 

9 

.31274  x  10-2 

.31284  x  lO2 

1  x  10-6 

10 

.28556  x  10-2 

.28574  x  10-2 

1.8  x  10-6 

As  we  can  see,  the  errors  in  the  k„'s  computed  using  the  numerical  method  are  slightly 
larger  than  the  previous  results  (Tables  1  through  4)  but  are  approximately  in  the  same 


order  of  magnitude. 
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Table  9.     POSITIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

Exact  k„  (m_1) 

WKB  Kn  (m-1) 

Absolute 
Error  (m-1) 

12 

.23190  x  10-2 

.23044  x  10-2 

1.5  x  10-5 

13 

.18524  x  ID"2 

.18494x  10-2 

3.0  x  10-6 

14 

.10756  x  10-2 

.10744  x  10-2 

1.2  x  10-6 

Table   10.     NEGATIVE  EXPONENTIAL  PROFILE,  TP  -  RB  MODES 


Mode 

Exact  k„  (m_1) 

WKB  Kn  (m-») 

Absolute 
Error  (m~l) 

23 

.82307  x  10~2 

.82309  x  10-2 

2  x  10-7 

24 

.79463  x  10-2 

.74469  x  10-2 

6x  10-7 

25 

.76664  x  10"2 

.76669  x  10-2 

5x  10-7 

26 

.73903  x  10~2 

.73909  x  10  2 

6  x  10-7 

27 

.71158  x  1<H 

.71 159  x  10-2 

1  x  10-' 

Table   11.     NEGATIVE  EXPONENTIAL  PROFILE,  PR  -  RB  MODES 


Mode 

Exact  xn  (nr1) 

WKB  K„  {nrl) 

Absolute 
Error  (nrl) 

30 

.61746  x  10-2 

.61689  x  10-2 

5.7  x  10-* 

31 

.57723  x  10-2 

.57689  x  10-2 

3.4  x  10-6 

32 

.530S2  x  10-2 

.53059  x  10-2 

2.3  x  10-6 

33 

.47670  x  10-2 

.47649  x  10-2 

2.1  x  10-6 
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IV.     CONCLUSIONS  AND  RECOMMENDATIONS 

The  WKB  approximations  of  acoustic  normal  modes  seem  to  give  results  accurate 
enough  to  be  of  practical  use.  Although  WKB  is  inherently  a  high  frequency  approxi- 
mation, meaning  that  it  works  better  for  higher  frequency,  the  results  from  our  tests 
show  that,  even  for  frequencies  around  one  Hertz,  this  technique  has  an  appreciable 
accuracy. 

When  applying  the  WKB  algorithm  for  arbitrary  sound  speed  profiles  the  determi- 
nation of  the  class  of  each  normal  mode  must  be  done  carefully.  Each  class  has  different 
formulae.  There  are  a  total  of  four  classes  in  a  single  duct  or  channel  environment. 

A  small  weakness  of  the  method  is  that  the  error  in  the  WKB  solution  increases 
when  a  turning  point  is  very  close  to  a  boundary.  However,  the  corresponding  errors  are 
still  very  small  for  the  exponential  and  hyperbolic  profiles  used  in  this  study. 

The  exact  analytical  solutions  to  the  Helmholtz  Equation  can  also  be  used  for 
comparison  to  any  other  methods.  These  exact  solutions  are  subject  to  pressure  release 
surface  and  rigid  bottom  boundary  conditions.  Exact  solutions  can  also  be  found  for 
any  boundary  conditions.  Specifically,  we  can  maintain  the  pressure  release  surface 
boundary  condition,  which  is  a  good  assumption,  and  use  a  non-rigid  bottom  boundary 
condition. 

Some  difficulties  were  found  when  working  with  Bessel  and  Associated  Legendre 
Functions.  The  IMSL  subroutines  used  for  Bessel  Function  evaluation  cannot  handle 
some  orders  and  values  of  the  argument.  With  respect  to  the  evaluation  of  Associated 
Legendre  Functions,  subroutines  are  not  available  in  the  IMSL  libraries.  Some  Fortran 
programs  were  coded  to  compute  them.  These  programs  also  work  only  for  certain 
ranges  of  order,  degree  and  arguments  of  the  functions.  Further  programming  work 
concerning  these  transcendental  functions  is  recommended. 

WKB  formulae  for  modes  trapped  in  the  water  column  over  a  non-rigid  bottom  were 
derived  in  Appendix  B  although  they  are  not  implemented  for  this  analysis.  A  general 
mathematical  expression  for  boundary  conditions  for  normal  modes  is  also  presented  in 
Appendix  B.  The  use  of  this  mathematical  expression  is  not  a  trivial  matter  because  the 
expression  depends  on  the  density  and  the  sound  speed  profile  of  the  sediment.  Further 
studies  on  using  WKB  algorithms  for  arbitrary  bottom  boundary  conditions  are  recom- 
mended. 
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In  this  thesis  only  single  duct/channel  environments  have  been  considered.  Essen- 
tially, we  have  ignored  double  duct/channel  problems.  However,  the  WKB  approxi- 
mation method  has  been  applied  to  other  nonacoustic  but  equivalent  double 
duct/channel  problems  (for  example,  transmission  of  electromagnetic  waves  through 
potential  barriers)  with  sucess  (Bender  and  Orszag,  1978). 

In  conclusion,  the  WKB  method  allows  for  accurate  and  fast  computation  of  normal 
modes  and  their  eigenvalues.  In  addition,  it  allows  for  storage  of  the  results  in  terms  of 
only  a  few  parameters  for  each  mode. 
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With 


APPENDIX  A.     FIRST  ORDER  \\  KB  THEORY 


d<f>  =  \x2n\  dz 


the  Helmholtz  Equation  becomes 


d2Z, 


i      dW\  ^l+Zh=s0  t 


d<p2       \k2„\      d($>      d(t> 


(A.l) 


This  equation  has  some  similarities  with  a  Bessel  Equation.  Let  us  try  a  solution  in  the 
form 


Zn  oc 


4> 


F(4>). 


Substituting  this  trial  solution  form  in  (A.l)  we  get 


dlF        1    dF 

dip2       <t>    d(\> 


i+r(z)- 


d/2f 
4>2     J 


F=0 


(A.2) 


with 


r{z)=\+ 


L  1  /    d\K:n\     \2_  J_  1  dl  I  Kzn  I 

4    IkJ4  V      dz      )"  2    |       |3"    ^ 


If  r(z)  is  negligible,  i.e.,  |r(z)|<^l,  (A.2)  can  be  approximated  as 


<TF        1     dF 

d<\>2      4>  d($> 


1- 


d/2r 

02 


F=0 


M-3) 


which  is  a  Bessel  Equation. 

The  general  solution  to  (A. 3),  for  k2„  >  0,  is 


F(<t>)  =  a'Jll2(<t>)+b'Yll2(<i>) 


{A  A) 


where  71/2  and  )'1/2  are  the  order  1/2  Bessel  Functions  of  the  l5'  and  2nd  kind,  respectively. 
But  since 
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,  sinx 

'1/2oc— — 


COS  X 

1/2  ~ 


(A. 4)  becomes 


and  thus 


V* 


sin  <b         cos  <b 
F(<}))  =  a—7±-+b—z+- 


a  sin  <f>+b  cos  </> 
Z»{z)  =  — 


K 


If  K\n  <  0,  instead  of  (A. 3)  we  would  get  the  Modified  Bessel  Equation  whose  general 
solution  is  given  by 

F(<f>)  =  a'Kll2(<l>)+b'Ill2{<f>) 

where  KU2  and  71/2  are  the  Modified  Bessel  Functions.   Since 


V 


0 


we  now  have 


71/2(0)  oc  — ==- 


_,  ,>       a?     +be 
F{4>)  = 


yj(t> 


and  thus 


VK 
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Let  us  assume  that  there  is  a  turning  point  at  z  =  0,  k22„>Q  for  z  >  0,  K*n<0  for 
2  <  0  and  near  the  turning  point  k2,„  =  yz.  So,  close  to  the  turning  point  we  have 


4>=\Z\K2n\dz  =  s\yXll\z\212 


with  s  =  zj\z\,  and  after  some  algebra, 

1  d  I  K:n 


Using  this  in  (A.l)  we  get 


K2n  |        <*/>  3</> 


</  Z„         i      </Z„ 


^02         30     i/0 
Let  us  now  look  for  a  solution  near  the  turning  point  in  the  form 


Substituting  this  solution  form  in  (A. 5)  we  get 


dlF        1     dF 


l- 


d/3r 

<p2 


F=0 


whose  general  solution  is 


F(0)  =  c15_1/3(c/))+C2J51/3(0) 


(A.5) 


where  B  represents  Bessel  Functions.  For  k2,„  >  0  they  are  7_]/3  and  J1/y  For  *£,  <  0  they 
are  7_1/3  and  71/3.  Relating  Z„  with  F  we  obtain 


zn  =  (-j<t>yi3^B_]l3(<j>)+c2Bll3m 


or 


v  1/3.     I 


ci#-i/ 


(±|v1/2UI3/2)+c2B1/3(±|y1/2UI3/2) 


(^■6) 
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where  the  plus  sign  is  for  k\„  <  0  and  the  minus  sign  for  K\n  >  0.  A  way  to  avoid  the  in- 
convenience of  sign  switching  is  to  use  Airy  Functions  which  are  related  to  B±hi  (see 
Chapter  III,  Section  A). 
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APPENDIX  B.     GENERAL  BOUNDARY  CONDITIONS 

A.     GENERAL  FORMULA 

Until  now  we  have  assumed  that  the  boundary  conditions  are  pressure  release  at  the 
surface  and  rigid  bottom.  The  first  one  is  a  good  assumption  but  the  second  one  can  be 
far  away  from  reality.  Let  us  now  present  a  general  bottom  condition  formula  for 
normal  modes  that  depends  on  the  bottom  characteristics. 

The  acoustic  pressure  due  to  the  n*  mode  is 

pn  =  Zn(z)Rn(r,e)elmt   .  (5.1) 

The  vertical  component  of  the  momentum  equation  relating  p„  to  the  vertical  particle 
velocity  uNn  in  the  n'h  mode  is 

where  p,  is  the  water  density  and  the  bottom  is  assumed  to  be  flat.   Using  (B.l)  in  (B.2) 
we  get 


p,a>      fc     Kn 


%„=„,„       ,-"V     •  (A3) 


At  the  bottom  boundary',  the  requirement  is  that  both  pn  and  uSn  are  continuous 
across  the  water-sediment  interface.  In  other  words,  the  normal  acoustic  impedance 
across  the  interface  is  continuous.  The  normal  acoustic  impedance  associated  with  the 
n'h  mode  is 

?«.-■&  ■  (5-4> 

By  using  (B.l)  and  (B.3)  in  (B.4)  one  can  obtain  as  a  general  condition 

/  dZn        p,co       \  ,  „  „% 

+/4L-Z„L//  =  0  (5.5) 


*  <*n 


if  zNn  at  the  water-sediment  interface  is  known  from  the  sediment  properties. 
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Let  us  use  a  local  plane  wave  approximation  of  the  normal  mode  near  z  =  II  in  the 
water  column,  i.e.. 

Zn(z)  =  An[e-iK»"{z-If)+ReiK»"(2-I[)] 

where  A„  is  a  constant,  R  is  the  complex  plane  wave  reflection  coefficient  and  KHn  is  the 
vertical  wavenumber  in  water  computed  at  the  boundary  (Clay  and  Medwin,  1977). 
With  this  plane  wave  approximate,  one  can  easily  show  that  (Kinsler  et  al.,  1982) 

P2c2b ,  D  ,v 

i      /     C'lb    W    C2b    \2    2 

where  the  index  1  refers  to  water,  2  refers  to  sediment  and  b  refers  to  the  value  at  the 
boundary.  R  is  also  known  as  the  Rayleigh  Reflection  Coefficient  which  can  be  equated 
as 


Pi  K2 


zn 


where  k2vi  is  the  vertical  wavenumber  in  the  sediment  layer  and  is  given  by 

2 

2         _w 2 

K2zn  ~       2        Kn     > 
C2 

and  Km  is  the  usual  vertical  wavenumber  in  the  water  column  (Kinsler  et  al.,  1982). 
Substituting  (B.6)  in  (B.5)  we  can  rewrite  (B.5)  as 

dZ„        p ,  \ 

-^+i-^K22nZnj2=H  =  0    .  (B.l) 

A  more  intuitive  form  for  (B.7)  is 

As  expected,  for  a  rigid  boundary  (R  =  1), 
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'dz"         =0 


dz 


r=H 


and,  for  a  pressure  release  boundary  (R  =  —  1), 

(z„W  =  o  . 

B.     TRAPPED  MODES  IN  THE  WATER  COLUMiN 

We  will  consider  only  modes  that  are  trapped  in  the  water  layer  in  this  WKB  for- 
mulation. This  means  that  K2zn  is  purely  imaginary  or 

K2zn  =  lfin 

with  /?„  real  and  defined  by 


/?„  =  -  K-4-    •  (*-9) 


Substituting  (B.9)  in  (B.7)  we  get 

(ifL-7r/?"z")=w=0  •  {DA0) 

Since  (B.7)  is  a  general  boundary  condition,  (B.10)  is  also  general  but  is  only  appli- 
cable to  modes  trapped  in  the  water  column.  We  will  apply  (B.10)  to  each  class  of 
normal  modes  in  the  water  column. 

1.     Class  I 

As  usual  we  assume  that  there  is  a  turning  point  at  z  =  z  with  0  <  z  <  H.  In  the 
exponential  region  (z  <  z  <  II)  the  solution  is 

Z„=     J cosh(0-fo)  (BM) 

V  I  Kzn  I 

with 


■=J> 


4>-\     \^zn\d2 


*2tl 


and 
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2     (0 

Kn  2 


As  (B.  1 1 )  must  satisfy  (B.  10),  4)b  must  be 


4>b  = 


H 


|  K2n  |  dz—  tanh" 


1         rfkrJ        Pi      A* 


2|ic„l 


2  dz 


p1      \K2n\ 


h=H 


(5.12) 


with  /?„  given  by  (B.9).  The  only  difference  relative  to  the  rigid  bottom  case  is  expression 
for  04.  All  the  other  formulae  remain  unchanged. 
2.     Class  II 

In  the  oscillatory  region,  z  <  z  <  H,  the  solution  is 


Z„  =  -=  sin(0+  -  )-  — —  cos(0+  —  j 

VK7* 


2  V 


(5.13) 


with 


0  =         I  K2  J 


</z 


and 


n 


4>s=         \Kzn\dz 


(B.13)  must  satisfy  (B.10),  implying  that 


sin<  (f)H+  -j-  +  tan  ' 


e*,+  »±Le-< 


e  *' 


-{D+Eje**  J 


=  0 


(5.14) 


where  D  is  defined  in  Chapter  II,  Equation  (22), 


<t> 


«mLK> 


,dz 


and 
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E  = 


Pi      Pn 

Pi     *zn    )Z=H 


(5.15) 


Therefore,  the  characteristic  equation  following  from  (B.14)  is 


lK^z==(n-i} 


for  D  <  0,  and 


I 


// 


Kzndz=  (n-—  )n-a 


for  D  >  0,  where 


a  =  tan 


<^r+ 


D+E    -<*    _, 


L  _l___(D+£)e^  J 


All  the  other  formulae  are  identical  to  the  rigid  bottom  case. 

3.  Class  III 

For  this  class  the  only  difference  relative  to  the  rigid  bottom  case  is  that  4>b  must 
now  be  computed  with  (B.12). 

4.  Class  IV 

The  solution  over  the  entire  ocean  depth  is 


1 


V  K2n 


sin  4> 


with 


dz 


To  satisfy  (B.10)  we  require 
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lK^Z  =  "n-iim~]{~DTE) 


for  D  >  0,  and 

■\(       1 


2ndz  =  nn+  tan 


for  D  <  0,  with  E  given  by  (B.15).  This  is  now  the  Class  IV  mode  characteristic  equation. 
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APPENDIX  C.     FORTRAN  PROGRAM  VVKBGEN 

A.  PROGRAM  DESCRIPTION 

The  program  WKBGEN  is  formed  by  a  main  program  and  several  subprograms. 
The  main  program  evaluates  the  minimum  sound  speed  and  then  controls  the  subpro- 
grams. The  inputs  must  be  given  in  the  PARAMETER  statement.  The  inputs  include 
the  ocean  depth  (H),  acoustic  frequency  (F),  the  wavenumber  increment  used  in  the  k„ 
search  (DK),  first  mode  (NMI)  and  last  one  (NMF)  to  be  computed.  The  units  in  this 
program  are  those  of  the  MKS  system.  The  input  sound  speed  profile  is  specified  using 
the  subprogram  FUNCTION  C(Z).  SUBROUTINE  TYPE  evaluates  the  class  of  each 
mode.  SUBROUTINE  CHARAC  controls  the  search  for  k„  and  the  turning  points.  The 
appropriate  characteristic  equation  is  given  by  FUNCTION  EQCHAR  and  the  turning 
points  are  evaluated  by  SUBROUTINE  ZTURN.  The  auxiliary  subprograms  FUNC- 
TION PHASE,  FUNCTION  FKZ2  and  FUNCTION  FKZ  evaluate  respectively  the 
phase  integral,  k\„  and  \kz„\.  The  Airy  Functions  are  computed  by  the  FUNCTION'S 
AI(Z)  and  BI(Z)  and  their  derivatives  by  FUNCTION'S  DA1(Z)  and  DBI(Z).  These  four 
subprograms  use  the  SUBROUTINE'S  MMBSJR  and  MMBSIR  in  the  IMSL  libraries 
to  evaluate  the  Bessel  and  Modified  Bessel  Functions,  respectively. 

As  output,  the  mode  numbers  and  the  respective  eigenvalues  are  printed  on  the 
screen.  These  results,  as  well  as  the  depths  of  the  turning  points,  are  also  written  in  a  file 
whose  name  is  specified  in  the  CALL  EXCMS  statement  at  the  begining  of  the  main 
program. 

The  numerical  method  used  to  solve  for  the  characteristic  equation  and  to  find  the 
turning  depths  is  the  simple  but  safe  Bisection  Method  (Gerald  and  Wheatley,  1989). 
The  method  to  compute  integrals  is  the  Trapezoidal  Rule  (Gerald  and  Wheatley,  1989). 
Derivatives  are  evaluated  by  forward  or  backward  finite  differences  (Gerald  and 
Wheatley,  1989). 

B.  PROGRAM  LISTING 

PROGRAM  WKBGEN 
C 
C 

c 

C     MAIN  PROGRAM 

C 

C 
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IMPLICIT  REAL*8  (A-H.O-Z) 

PARAMETER  (H=        ,F=      ,DK=      ,NMI=   ,NMF=   ) 


CALL  EXCMSC'FILEDEF  1  DISK  WKBGEN  DATA  A') 

OM=8.  D0*DATAN( 1.  D0)*F 

ZEX=DMOD(H,10.DO) 


CS=C(0) 
CB=C(H) 


CM=99999.DO 

DO  Z=O.DO,H-ZEX,10.DO 

CM=DMIN1(CM,C(Z)) 
END  DO 
CM=DMIN1(CM,C(H)) 


XK0=OM/CM 
XKL=XK0 


DO  N=NMI,NMF 

DO  XKI=XKL,0,DK 

XKF=XKI+DK 
CNI=OM/XKI 
CNF=OM/XKF 

CALL  TYPE(CNI,CS,CB,NTYPE) 
CALL  ZTURN(CNI,H,ZEX,NTYPE,ZI1,ZI2) 
CALL  ZTURN(CNF,H,ZEX,NTYPE,ZF1,ZF2) 

CALL  CHARAC(OM,H,ZEX,N,XKI,XKF,NTYPE,ZIl,ZI2,ZFl,ZF2,NSOL,XKSOL, 
$  ZT1,ZT2) 

IF  (NSOL.  EQ.  1)  GO  TO  100 

END  DO 

100   PRINT*, N,XKSOL 

IF  (NTYPE.EQ.  l.OR.  NTYPE.EQ.  2)  THEN 

WRITE  (1,*)  N,NTYPE,XKS0L,ZT1 
ELSE  IF  (NTYPE.EQ. 3)  THEN 

WRITE  (1,*)  N,NTYPE,XKS0L,ZT1,ZT2 
ELSE 

WRITE  (1,*)  N,NTYPE,XKSOL 
END  IF 

XKL=XKSOL 
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END  DO 

STOP 

END 
C 
C 

c 

C 

c 
c 

SUBROUTINE  TYPE(CN,CS,CB,NTYPE) 
REAL*8  CN,CS,CB 
C 

IF  (CN.  GT.CS.AND.CN.  LT.  CB)  THEN 

NTYPE=1 
ELSE  IF  (CN.  LT.CS.AND.CN.  GT.  CB)  THEN 

NTYPE=2 
ELSE  IF  (CN.  LT.CS.AND.CN.  LT.  CB)  THEN 

NTYPE=3 
ELSE 

NTYPE=4 
END  IF 
C 

RETURN 
END 
C 
C 

c 

C 
C 
C 

SUBROUTINE  CHARAC(OM>H,ZEX,NM,XI,XFJNTJZIl,ZI2,ZFl,ZF2,NSOL,XSOL) 
$  Z21.Z22) 

IMPLICIT  REAL*8  (A-H,0-Z) 
XSOL=9999.D0 
NSOL=l 

FI=EQCHAR(NM,NT,H,OMJXI,ZIl,ZI2) 
FF=EQCHAR( NM , NT , H , OM , XF , ZF 1 , ZF2 ) 
IF  ((FI*FF).GT.  0.D0)  THEN 

NSOL=0 

GO  TO  100 
ELSE  IF  (FI.EQ.  0.  DO)  THEN 

XSOL=XI 

GO  TO  500 
ELSE  IF  (FF.  EQ.  0.D0)  THEN 

XSOL=XF 

GO  TO  500 
ELSE 

CONTINUE 
END  IF 

DO  200  N=l,500 
XM2=(XF+XI)/2.D0 
IF(N.EQ.  1)  GO  TO  50 
IF  ((XM2-XM1).EQ.  0.D0)  THEN 

XSOL=XM2 
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GO  TO  500 
END  IF 
50  XM1=XM2 

CN2=OM/XM2 

CALL  ZTURN(CN2,H,ZEX,NT,Z21,Z22) 

F2=EQCHAR(NM,NT,H,OM,X2,Z21,Z22) 

IF  ((F2*FF). LT. 0.  DO)  THEN 

XI=XM2 
ELSE 

XF=XM2 
FF=F2 
END  IF 
200  CONTINUE 
100  CONTINUE 
500  CONTINUE 
RETURN 
END 
C 
C 

c 

c 
c 
c 

SUBROUTINE  ZTURN(CN,H, ZEX.NT, ZT1 ,ZT2) 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*8  NM 

ZT1=99999.D0 

ZT2=99999.D0 

DO  100  D=  0.D0,H,10.D0 

XF=D+10.  DO 

XI=D 

FI=C(XI)-CN 

FF=C(XF)-CN 

IF  ((FI*FF).GT.  0.D0)  THEN 

GO  TO  100 
ELSE  IF  (FI.EQ.  O.DO)  THEN 

XSOL=XI 

GO  TO  500 
ELSE  IF  (FF.EQ.  O.DO)  THEN 

XSOL=XF 

GO  TO  500 
ELSE 

CONTINUE 
END  IF 

DO  200  N=l,500 
XM2=(XF+XI)/2.D0 
IF(N.EQ.  1)  GO  TO  50 
IF  ((XM2-XM1).EQ.  O.DO)  THEN 

XS0L=XM2 

GO  TO  500 
END  IF 
50  XM1=XM2 

F2=C(XM2)-CN 

IF  ((F2*FF).LT.  O.DO)  THEN 

XI=XM2 
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ELSE 

XF=XM2 

FF=F2 

END  IF 

200  CONTINUE 

100  CONTINUE 

500  CONTINUE 

ZT1=XS0L 

IF  (NT.  EQ.  3)  THEN 
DO  110  D=  H,0.D0,-10.D0 
XF=D-10.D0 
XI=D 

FI=C(XI)-CN 
FF=C(XF)-CN 
IF  ((FI*FF).GT.  0.D0)  THEN 

GO  TO  110 
ELSE  IF  (FI.EQ.  0.D0)  THEN 

XSOL=XI 

GO  TO  510 
ELSE  IF  (FF.EQ.  0.D0)  THEN 

XSOL=XF 

GO  TO  510 
ELSE 

CONTINUE 
END  IF 

DO  210  N=l,510 
XM2=(XF+XI)/2.D0 
IF(N.  EQ.  1)  GO  TO  51 
IF  ((XM2-XM1).EQ.  O.DO)  THEN 

XS0L=XM2 

GO  TO  510 
END  IF 
51  XM1=XM2 

F2=C(XM2)-CN 

IF  ((F2*FF).LT.  O.DO)  THEN 

XI=XM2 
ELSE 

XF=XM2 

FF=F2 
END  IF 
210  CONTINUE 
110  CONTINUE 
510  CONTINUE 
ZT2=XS0L 
END  IF 
RETURN 
END 
C 
C 

c 

C 
C 

c 

FUNCTION  EQCHAR(NM,NT,H,0M,XI,ZT1,ZT2) 
IMPLICIT  REAL*8  (A-H,0-Z) 


58 


PI=4.D0*DATAN(1.D0) 

D=(FKZ(OM,XI,H)-FKZ(OM,XI,H-5.D0))/(10.D0*FKZ(OM,XI,H)**2) 
IF  (DABS(D).GT.  .  999D0)  THEN 

D1=(D/DABS(D))*.  999D0 
ELSE 

D1=D 
END  IF 


IF  (NT.  EQ.  1)  THEN 

EQCHAR=PHASE(0.D0,ZTl,OM,XI)-(NM-.25D0)*PI 
IF  (H-ZT1.GT.  10.D0)  THEN 

PHIB=PHASE(ZT1,H,0M,XI)-DATANH(D1) 

EQCHAR=EQCHAR+DATAN(DEXP( -2.  D0*PHIB)/2.  DO) 
ELSE 

GAM=(FKZ2(0M,XI,H)-FKZ2(0M,XI,ZT1))/(ZT1-H) 

ZH=(GAM/DABS(GAM))*DABS(GAM)**(1.D0/3.D0)*(H-ZT1) 

A=DBI(ZH) 

B=DAI(ZH) 

EQCHAR=EQCHAR-DATAN( B/A) 
END  IF 


ELSE  IF  (NT. EQ. 2)  THEN 

EQCHAR=PHASE ( ZT1 , H , OM , XI ) 
IF(ZT1.GT.  10.  DO)  THEN 
PHIS=PHASE( 0.  DO , ZT1 , OM , XI ) 
ALF1=DEXP(PHIS)+D*DEXP(-1.D0*PHIS)/2.D0 
ALF2=((DEXP(-1.D0*PHIS)/2.D0)-D*DEXP(PHIS)) 
IF  (DLOG10(DABS(ALF1))-DLOG10(DABS(ALF2)).GT.  60.  DO)   THEN 

IF  (D. GE. 0)  THEN 
ALF=PI/2.  DO 

ELSE 

ALF=-1.D0*PI/2.D0 

END  IF 
ELSE 

ALF=DATAN ( ALF 1 / ALF  2 ) 
END  IF 
IF  (D.GE.O.DO)  THEN 

EQCHAR=EQCHAR-(NM-1.  25D0)*PI+ALF 
ELSE 

EQCHAR=EQCHAR-(NM-.  25D0)*PI+ALF 
END  IF 
ELSE 
•   GAM=(FKZ2(OM,XI,ZT1)-FKZ2(OM,XI,O.DO))/ZT1 
ZS=(GAM/DABS(GAM) )*DABS(GAM)**( 1.  DO/3.  DO) 
A=BI(ZS) 
B=AI(ZS) 
IF  (D.LT.  O.DO)  THEN 

EQCHAR=EQCHAR-(NM-.25D0)*PI+DATAN((A+D*B)/(B-D*A)) 
ELSE 

EQCHAR=EQCHAR-(NM-1.25D0)*PI+DATAN((A+D*B)/(B-D*A)) 
END  IF 
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END  IF 


ELSE  TF  (NT.  EQ.  3)  THEN 
IF  (H-ZT2.GT.  10.  DO)  THEN 

PHIB=PHASE( ZT2 , H , OM , XI ) -DATANH( Dl ) 

ALF2=DATAN(DEXP( -2.  D0*PHIB)/2.  DO) 
ELSE 

GAM=(FKZ2(OM,XI,H)-FKZ2(OM,XI,ZT2))/(ZT2-H) 

ZH=(GAM/DABS(GAM))*DABS(GAM)**(1.D0/3.D0)*(H-ZT2) 

A=DBI(ZH) 

B=-l.DO*DAI(ZH) 

ALF2=DATAN(B/A) 
END  IF 
IF(ZT1.GT.  10.  DO)  THEN 

PHIS=PHASE(O.DO,ZT1,OM,XI) 

ALF1=DATAN(DEXP( -2.  D0*PHIS)/2.  DO) 
ELSE 

GAM=( FKZ2( OM , XI , ZT1 ) -FKZ2( OM , XI , 0.  DO) ) /ZT1 

ZS=(GAM/DABS(GAM))*DABS(GAM)**(1.D0/3.D0)*ZT1 

A=BI(ZS) 

B=AI(ZS) 

ALF1=DATAN(B/A) 
END  IF 

EQCHAR=PHASE( ZT1 , ZT2 , OM , XI ) 

EQCHAR=EQCHAR-(NM-.5D0)*PI-ALF1+ALF2 


ELSE 

EQCHAR=PHASE(0.D0,H,OM,XI) 

IF  (DABS(D).GT.  l.D-60)  THEN 

E=1.D0/D 

IF(D.GT.  O.DO)  THEN 

EQCHAR=EQCHAR-NM*PI+DATAN( E ) 
ELSE 

EQCHAR=EQCHAR-NM*PI-DATAN(E) 
END  IF 
ELSE 

EQCHAR=EQCHAR-NM*PI+PI/2.  DO 
END  IF 


END  IF 


RETURN 

END 
C 
C 

c 

c 
c 
c 
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FUNCTION  PHASE(A,B,OM,XKN) 

IMPLICIT  REAL*8  (A-H,0-Z) 

XNL=DINT((B-A)/5.D0) 

NL=NINT(XNL) 

EX=DM0D((B-A),5.D0) 

PHASE=0.  DO 

DO  1=1, NL 

XI=DBLE(I) 

PHASE=PHASE+( FKZ( OM , XKN , ( XI - 1.  DO )*5 .  DO+A) 
$         +FKZ(OM, XKN, XI*5.  DO+A) )*2.5D0 
END  DO 

PHASE=PHASE+( FKZ ( OM , XKN , XNL*5 .  DO+A) 
$      +FKZ(0M,XKN,B))*EX/2.D0 
RETURN 
END 
C 
C 

c 

C 

c 
c 

FUNCTION  FKZ(OM,XKN,Z) 

REAL* 8  FKZ,OM,XKN,Z 

FKZ=(0M/C(Z))**2-XKN**2 

FKZ=DSQRT( DABS( FKZ ) ) 

RETURN 

END 
C 
C 

c 

c 
c 
c 

FUNCTION  FKZ2(OM,XKN,Z) 

REAL* 8  FKZ2,OM,XKN,Z 

FKZ2=(OM/C(Z))**2-XKN**2 

RETURN 

END 
C 
C 

c 

C 
C 
C 

FUNCTION  C(Z) 
REAL*8  C,Z 
C= 

RETURN 
END 
C 

c 
c 
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FUNCTION  AI(Z) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J13,JM13,I13,IM13 
DIMENSION  RJ(2),WK(4),B(2) 


IF  (Z.EQ.  0.D0)  THEN 

AI=. 35502805D0 

RETURN 
ELSE  IF  (Z.LT.  O.DO)  THEN 


ARG=2. DO*(DSQRT(  -1.  D0*Z)**3)/3.  DO 
N=2 


ORDER=2.D0/3.D0 

CALL  MMBS JR(ARG, ORDER, N,RJ,WK,IER) 

JM13=4.D0*RJ(1)/(3.D0*ARG)-RJ(2) 


ORDER=l.  DO/3.  DO 

CALL  MMBS JR(ARG, ORDER, N,RJ,WK,IER) 
J13=RJ(1) 


AI=DSQRT( -1. DO*Z)*( JM13+J13)/3.  DO 


RETURN 
ELSE 


ARG=2.  D0*(DSQRT(Z)**3)/3.  DO 

0RDER=2.D0/3.D0 

NB=2 

IOPT=l 


CALL  MMBSIR(ARG,ORDER,NB,IOPT,B,IER) 
IM13=4. DO*B( l)/(3. D0*ARG)+B(2) 


ORDER=l.D0/3.D0 
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CALL  MMBSIR(ARG, ORDER, NB , IOPT,B , IER) 
I13=B(1) 


AI=DSQRT(Z)*(IM13-I13)/3.D0 
RETURN 


END  IF 


C 
C 
C 

C 
C 
C 

END 
C 
C 
C 

C 

c 
c 

FUNCTION  BI(Z) 
IMPLICIT  REAL*8  (A-H.O-Z) 
REAL*8  J13,JM13,I13,IM13 
DIMENSION  RJ(2),WK(4),B(2) 

C 

C 

C 


IF  (Z.EQ.  0.D0)  THEN 

BI=. 61492663D0 

RETURN 
ELSE  IF  (Z.  LT.  O.DO)  THEN 


ARG=2. DO*(DSQRT( -1.  D0*Z)**3)/3.  DO 

N=2 


ORDER=2.D0/3.D0 

CALL  MMBS JR(ARG, ORDER, N,RJ,WK, IER) 

JM13=4.D0*RJ(1)/(3.D0*ARG)-RJ(2) 


ORDER=l.D0/3.D0 

CALL  MMBS JR(ARG, ORDER, N,RJ,WK, IER) 

J13=RJ(1) 


BI=DSQRT( -1. D0*Z/3.  D0)*( JM13-J13) 
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RETURN 
ELSE 


ARG=2.D0*(DSQRT(Z)**3)/3.D0 

ORDER=2.D0/3.D0 

NB=2 

I0PT=1 


CALL  MMBSIR(ARG, ORDER, NB , IOPT,B , IER) 
IM13=4.  DO*B(  l)/(3.  D0*ARG)+B(2) 


ORDER=l.D0/3.D0 

CALL  MMBSIR(ARG, ORDER, NB,IOPT,B, IER) 

I13=B(1) 


BI=DSQRT(Z/3.D0)*(IM13+I13) 
RETURN 


END  IF 


C 
C 

c 

c 
c 
c 

END 
C 

c 
c 

c 
c 
c 

FUNCTION  DAI(Z) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J23,JM23,I23,IM23 
DIMENSION  RJ(2),WK(4),B(2) 

C 

C 

C 


IF  (Z.EQ.  O.DO)  THEN 
DAI=-. 25881940D0 
RETURN 

ELSE  IF  (Z.  LT.  O.DO)  THEN 
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ARG=2.D0*(DSQRT(-1.  D0*Z)**3)/3.  DO 
N=2 


ORDER=l.  DO/3.  DO 

CALL  MMBS JR( ARG , ORDER , N , RJ , WK , IER) 
JM23=2.D0*RJ(1)/(3.D0*ARG)-RJ(2) 


ORDER=2.D0/3.DO 

CALL  MMBS JR( ARG, ORDER, N,RJ,WK, IER) 

J23=RJ(1) 


DAI=Z*(JM23-J23)/3.  DO 


RETURN 
ELSE 


ARG=2. DO*(DSQRT(Z)**3)/3.  DO 
ORDER=l.D0/3.D0 

NB=2 
IOPT=l 


CALL  MMBSIR(ARG,ORDER,NB,IOPT,B,IER) 
IM23=2.D0*B(1)/(3.D0*ARG)+B(2) 


ORDER=2.D0/3.D0 

CALL  MMBSIR( ARG, ORDER, NB,IOPT,B, IER) 

I23=B(1) 


DAI=-1. D0*Z*(IM23-I23)/3.  DO 
RETURN 


END  IF 


END 
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c 

c 
c 
c 

FUNCTION  DBI(Z) 
IMPLICIT  REAL*8  (A-H,0-Z) 
REAL*8  J23,JM23,I23,IM23 
DIMENSION  RJ(2),WK(4),B(2) 

C 

C 

c 


IF  (Z.EQ.  O.DO)  THEN 

DBI=.  44828836D0 

RETURN 
ELSE  IF  (Z.LT.  O.DO)  THEN 


ARG=2. DO*(DSQRT( -1.  D0*Z)**3)/3.  DO 
N=2 


ORDER=l.  DO/3.  DO 

CALL  MMBS JR(ARG, ORDER, N,RJ,WK,IER) 
JM23=2. DO*RJ( 1 ) /( 3.  DO*ARG) -RJ( 2 ) 


ORDER=2.D0/3.D0 

CALL  MMBSJR(ARG,ORDER,N,RJ,WK,IER) 

J23=RJ(1) 


DBI=-1. DO*Z*( JM23+J23)/DSQRT(3.  DO) 


RETURN 
ELSE 


ARG=2. D0*(DSQRT(Z)**3)/3.  DO 

0RDER=1.D0/3.D0 

NB=2 

IOPT=l 


CALL  MMBSIRCARG, ORDER, NB,IOPT,B,IER) 

IM23=2.D0*B(1)/(3.D0*ARG)+B(2) 
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ORDER=2.D0/3.D0 

CALL  MMBS IR( ARG , ORDER , NB , IOPT , B , IER) 

I23=B(1) 


DBI=Z*( IM23+I23)/DSQRT( 3.  DO) 
RETURN 


END  IF 


END 
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P54262  Pimentel 

c.l       Computation  of  acous- 
tic normal  modes  in  the 
ocean  using  asymptotic 
expansion  methods. 


